- 77n- iOD-, 730 


NASA Technical Memorandum 1 02730 


NASA-TM- 102730 19910002360 


DEVELOPMENT OF UNSTRUCTURED GRID METHODS FOR STEADY 
AND UNSTEADY AERODYNAMIC ANALYSIS 


JOHN T. BATINA 


NOVEMBER 1990 


fW\SA 

National Aeronautics and 
Space Administration 

Langley Research Center 

Hampton, Virginia 23665 


o-'.r, 

fc S v tT. / ♦ ■ • / r ■ *■=•?-' *, 

l-tk i e -it's; „ (■ r < - 

::ov o • 1990 



ii 


U 13 RAW NASA 

HA’.^ION, vikc:,ni 




N 



DEVELOPMENT OF UNSTRUCTURED GRID METHODS FOR 
STEADY AND UNSTEADY AERODYNAMIC ANALYSIS 


John T. Batina* 

NASA Langley Research Center 
Hampton, Virginia 23665-5225 USA 


Abstract 

The current status of the development of unstructured grid 
methods in the Unsteady Aerodynamics Branch at NASA Langley 
Research Center is described. These methods are being 
developed for steady and unsteady aerodynamic applications. 
The paper first highlights the flow solvers that have been 
developed for the solution of the unsteady Euler and Navier- 
Stokes equations and then gives selected results which 
demonstrate various features of the capability. The results 
demonstrate two- and three-dimensional applications for both 
steady and unsteady flows. Comparisons are also made with 
solutions obtained using a structured grid code and with 
experimental data to determine the accuracy of the 
unstructured grid methodology. These comparisons show good 
agreement which thus verifies the accuracy. 


Introduction 

Considerable progress has been made over the past two 
decades on developing computational fluid dynamics (CFD) 
methods for aerodynamic analysis. “1 .2 Recent work in CFD has 
focused primarily on developing algorithms for the solution of 
the Euler and Navier-Stokes equations. For unsteady 
aerodynamic and aeroelastic analysis, these methods generally 
require that the mesh move to conform to the instantaneous 
position of the moving or deforming body under consideration. 
Many of the methods that are currently being developed assume 
that the mesh moves rigidly or that the mesh shears as the 
body deforms. These assumptions consequently limit the 
applicability of the procedures to rigid-body motions or 
small-amplitude deformations. Furthermore, these methods 
of solution typically assume that the computational grid has an 
underlying geometrical structure. As an alternative, 
algorithms have been developed recently which make use of 
unstructured grids.3-16 j n two dimensions these grids are 
typically made up of triangles and in three dimensions they 
consist of an assemblage of tetrahedra. The unstructured grid 
methods have several distinct advantages over structured grid 
methods which make them attractive for unsteady aerodynamic 
and aeroelastic analyses. For example, the primary 
advantage of the unstructured grid methodology is the ability to 
easily model very complicated three-dimensional 
configurations with virtually unlimited geometrical 
complexity. A second advantage is that the methodology allows 
for a general way to move the mesh to treat realistic motions 
and structural deformations of complete aircraft 
configurations. The deforming grid capability does not involve 
any assumptions which limit applications to small 
deformations, such as simple grid shearings which are done in 
some structured grid codes. A third advantage is that it enables 
in a natural way for adaptive mesh refinement to more 
accurately predict the physics of the flow. Highly accurate 
solutions may thus be obtained by using far fewer grid points 
than if a globally fine mesh was used. 

The purpose of the paper is to describe the current status 
of the development of unstructured grid methods within the 
Unsteady Aerodynamics Branch at NASA Langley Research 
Center. 10-16 jh 0 paper first highlights the flow solvers that 
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have been developed for solution of the time-dependent Euler 
and Navier-Stokes equations and then gives selected results 
which demonstrate various features of the capability. The flow 
solvers that are described are either of the central- 
difference-type with explicit artificial dissipation or of the 
upwind-type which are naturally dissipative. Both implicit 
and explicit temporal discretizations are discussed for the 
time-integration of the governing fluid flow equations. Details 
on the adaptive mesh refinement procedures are also given. 
The selected results that are presented demonstrate two- and 
three-dimensional applications for both steady and unsteady 
flows. Comparisons are also made with solutions obtained 
using a structured grid code and with experimental data to 
determine the accuracy of the unstructured grid methodology. 


Central- Diff erence-Type _Flow Solver 

The unsteady Euler equations in integral form are solved 
using a finite-volume algorithm that was developed for use on 
unstructured grids of triangles in 2D or tetrahedra in 
3 D J 0-1 3 j he algorithm reduces conceptually to central 
differencing on a rectangular mesh and thus is referred to as a 
central-difference-type flow solver. With this solver, 
artificial dissipation is added explicitly to prevent oscillations 
near shock waves and to damp high-frequency uncoupled error 
modes. Specifically, an adaptive blend of harmonic and 
biharmonic operators is used, corresponding to- second and 
fourth difference dissipation, respectively. The biharmonic 
operator provides a background dissipation to damp high 
frequency errors and the harmonic operator prevents 
oscillations near shock waves. 

The Euler equations are integrated in time using a 
standard, explicit, four-stage, Runge-Kutta time-stepping 
scheme. In this scheme the convective operator is evaluated at 
each stage and, for computational efficiency, the dissipative 
operator is evaluated only at the first stage. The scheme is 
second-order-accurate in time and includes the necessary 
terms to account for changes in cell volumes due to a moving or 
deforming mesh. Furthermore, this explicit scheme has a step 
size that is limited by the Courant-Friedricks-Lewy (CFL) 

condition corresponding to a CFL number of 2 ^ 2 . To accelerate 
convergence to steady-state, the CFL number may be increased 
by averaging implicitly the residual with values at 
neighboring grid points. These implicit equations are solved 
approximately using several Jacobi iterations. Convergence to 
steady-state is further accelerated using enthalpy damping and 
local time stepping. The local time stepping uses a maximum 
allowable step size at each grid point as determined by a local 
stability analysis. For unsteady applications, however, a 
global time step is usually used because of the time-accuracy 
requirement. The maximum allowable global time step may be 
increased to a value that is larger than that dictated by the CFL 
condition, by using a time accurate version of the residual 
smoothing. 

If interest is restricted to supersonic flow past conical 
bodies, then the conical flow assumption can be made. This 
reduces the problem from three dimensions to two dimensions, 
which significantly decreases the computational resources that 
are required to investigate such flows. The conical flow 
assumption is exact for inviscid supersonic flow and is used to 
efficiently investigate vortex-dominated flows. For viscous 





flow, however, a length dependence remains in the Reynolds 
number (Re), although the flow may be considered to be 
locally conical. The Reynolds number therefore determines 
the location of the plane at which the solution is determined. 
Conical Euler and Navier-Stokes solvers have been developed 
based on a central-difference-type spatial discretization, 1 3 
similar to the two- and three-dimensional Euler solvers 
described above. The viscous fluxes in the conical Navier- 
Stokes solver are evaluated by first computing derivatives of 
the velocity components required by the shear stresses and 
heat flux terms using Green's theorem. Once evaluated, these 
fluxes are averaged across edges in the same manner as the 
inviscid fluxes. Also, the viscous fluxes are computed only at 
the first stage of the Runge-Kutta time-marching scheme for 
computational efficiency. This effectively reduces the 
computational work to evaluate the viscous fluxes by a factor 
of four. 


A refinement indicator is used to determine if an element 
in the mesh is to be refined or subdivided into smaller 
elements. Typically, the absolute change in density along an 
edge is used as an indicator for flows with shock waves and 
total pressure loss is used for flows with vortices. The 
refinement indicator is compared with a preset tolerance to 
determine whether a given element should be refined. If the 
tolerance is exceeded, a new node is created at the midpoint of 
the edge and the element is divided. Each time the mesh is 
refined, an element may be divided in one of several different 
ways. The coordinates of the new node are determined by 
averaging the coordinates of the endpoints that make up tho 
bisected edge. Special care must be taken, however, when an 
edge that is to be divided lies on a boundary of the grid, since 
the midpoint of the edge does not generally lie on the boundary. 
In this case, the location of the new node is determined 
generally by using a spline of the boundary coordinates. 


Upwind-Type Flow. Solver 


Def orming Me s h Algorithm 


The unsteady Euler equations may be solved alternatively 
by using upwind differencing and either flux-vector or flux- 
difference splitting similar to upwind schemes developed for 
use on structured meshes. 1 °» 1 4-1 6 The present 
unstructured grid algorithm is thus referred to as an upwind- 
type flow solver. The spatial discretization of this solver 
involves a so-called flux-split approach based on either the 
flux-vector splitting of van Leer 17 or the flux-difference 
splitting of Roe. 18 These flux-split discretizations account 
for the local wave-propagation characteristics of the flow and 
they capture shock waves sharply with at most one grid point 
within the shock structure. A further advantage is that these 
discretizations are naturally dissipative and consequently do 
not require additional artificial dissipation terms or the 
adjustment of free parameters to control the dissipation. 
However, in calculations involving higher-order upwind 
schemes such as these, oscillations in the solution near shock 
waves are expected to occur. To eliminate these oscillations 
flux limiting is usually required. In the present study, a 
continuously differentiable flux limiter was employed. 14 ” 18 

The Euler equations are integrated in time using either an 
explicit Runge-Kutta method (described in the previous 
section) or an implicit time-integration scheme involving a 
Gauss-Seidel relaxation procedure. 14 The procedure is 
implemented by re-ordering the elements that make up the 
unstructured mesh from upstream to downstream. The 
solution is obtained by sweeping two times through the mesh as 
dictated by stability considerations. The first sweep is 
performed in the direction from upsteam to downstream and 
the second sweep is from downstream to upstream. For purely 
supersonic flows the second sweep is unnecessary. This 
relaxation scheme is stable for large time steps and thus 
allows the selection of the step size based on the temporal 
accuracy of the problem being considered, rather than on the 
numerical stability of the algorithm. Consequently, very large 
time steps may be used for rapid convergence to steady state, 
and an appropriate step size may be selected for unsteady 
cases, independent of numerical stability issues. 


Spatial Adaption . Procedure 

Spatial adaption is employed with the unstructured grid 
flow solvers to enrich the mesh locally in regions of high 
spatial flow gradients to more accurately and efficiently 
resolve the physics of the flow. 13 Equally attractive are 
coarsening techniques that remove elements from regions 
where relatively small changes in the flow variables occur. 
Both enrichment and coarsening procedures are currently 
being developed. However, only the enrichment procedure is 
described as follows. 


For problems where the aircraft moves or deforms, the 
mesh must move so that it continuously conforms to the 
instantaneous shape or position of the vehicle. This is 
accomplished by using a spring network to model the original 
mesh such that each edge of the triangle or tetrahedron is 
represented by a sprihg. 12 The spring stiffness for a given 
edge is taken to be inversely proportional to the length of the 
edge. Grid points on the outer boundary of the mesh are held 
fixed and the instantaneous locations of the points on the inner 
boundary (aircraft) are given by the prescribed surface 
motion. At each time step, the static equilibrium equations in 
the x, y, and z directions, which result from a summation of 
forces, are solved iteratively at each interior node of the grid 
for the displacements. This is accomplished by using a 
predictor-corrector procedure, which first predicts the 
displacements of the nodes by extrapolation from grids at 
previous time levels and then corrects these displacements 
using several Jacobi iterations of the static equilibrium 
equations. The predictor-corrector procedure has been found 
to be more efficient than simply performing Jacobi iterations 
because far fewer iterations are required to achieve acceptable 
convergence. In practice, it has been found that only one or 
two iterations are sufficient to accurately move the mesh. 


ResultS-and-Discussica 

Selected results from the unstructured grid methods of 
Refs. 10-16 are presented for two- and three-dimensional 



Fig. 1 Partial view of unstructured grid of triangles about 
the NACA 0012 airfoil. 
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geometries for both steady and unsteady flows. Comparisons 
are made with solutions obtained using a structured grid code 
and with experimental data to determine the accuracy of the 
methodology. 

Two-Dimens ional Euler Results 

To assess the two dimensional central-difference and 
upwind-type Euler solvers, calculations were performed for 
the NACA 0012 airfoil. These results were obtained by using 
the unstructured grid shown in Fig. 1 which was generated 
using the advancing front method. 6 * 7 The grid has 3300 
nodes, 6466 triangles, and extends 20 chordlengths from the 
airfoil with a circular outer boundary. Also there are 110 
points that lie on the airfoil surface. This is the same mesh 
that was used to obtain the results that were presented in Refs. 
10, 11, and 14. 

Generalized aerodynamic forces were obtained using the 
central-difference-type Euler solver for the NACA 0012 
airfoil oscillating in either plunge or pitch-about-the- 
quarter-chord to determine the accuracy of the unstructured 
grid method for unsteady aerodynamic applications. 
Calculations were performed for the airfoil at a freestream 


Mach number of » 0.8 and zero degrees angle of attack. 
Comparisons are made with parallel computations performed 
using the CFL3D 1 9 code, run in a 2-D mode. The CFL3D code 
is an Euler/Navier-Stokes code based on cell-centered 
upwind-difference discretizations of the governing flow 
equations based on structured meshes. The CFL3D Euler 
results that were used in the present study were obtained from 
Ref. 20. Generalized aerodynamic forces for the NACA 0012 
airfoil are presented in Fig. 2, computed using the so-called 
pulse analysis within the unstructured grid Euler code, the 
unstructured grid Euler harmonic analysis, and the CFL3D 
harmonic analysis. The results are plotted as real and 

imaginary components of the unsteady forces, A Uj , as a 

function of reduced frequency k. Plunge and pitch-about-the- 
quarter-chord motions are defined as modes 1 and 2, 
respectively. Thus, for example, A 12 is the lift coefficient due 
to pitching. Both sets of harmonic results were obtained at six 
values of reduced frequency: k = 0.0, 0.125, 0.25, 0.5, 0.75, 
and 1.0. As shown in Fig. 2, the forces from the unstructured 
grid Euler pulse analysis agree well with the forces from the 
harmonic analysis. The harmonic analysis, however, is 
considered to be the more accurate of the two sets of 
calculations, since the local linearity assumption in the pulse 
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Fig. 2 Comparisons of generalized aerodynamic forces 

for the NACA 0012 airfoil at M_= 0.8 and 0° 
computed using CFL3D and the unstructured-grid 
central-difference-type Euler flow solver. 


Fig. 3. Comparisons of steady-state results for the NACA 
0012 airfoil at M_= 0.8 and «„» 1.25° computed 
using the upwind-type Euler flow solver with flux- 
vector splitting. 
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analysis may be questionable for transonic flow cases. As 
further shown in Fig. 2, the harmonic forces from the 
unstructured grid Euler code are in very good agreement with 
the CFL3D (Euler) harmonic forces, for both plunge and pitch 
motions, for all values of reduced frequency that were 
considered. 

To test the more-recently-developed upwind-type Euler 
flow solver, steady flow results were obtained for the NACA 

0012 airfoil at = 0.8 and a o = 1.25°, using both 

implicit and explicit time-marching. The explicit time- 
marching results were obtained using a CFL number of 2.5 
(since the CFL limit is approximately 2.8) and the implicit 
time-marching results were obtained using a CFL number of 
100,000. Such a large value was used for the implicit results 
since the relaxation scheme has maximum damping and hence 
fastest convergence for very large time steps. This is in 
contrast with implicit approximate factorization schemes 
which have maximum damping for CFL numbers on the order 
of 10. A comparison of the convergence histories between 
explicit and implicit time-marching is shown in Fig. 3(a). 
The "error” in the solution was taken to be the L 2 -norm of the 
density residual. As shown in Fig. 3(a), the explicit solution 
is very slow to converge. This solution takes approximately 
10,000 time steps to become converged to engineering 
accuracy, which is taken to be a four order of magnitude 
reduction in solution error. In contrast, the implicit solution 
is converged to four orders of magnitude in only approximately 
500 steps and is converged to machine zero is less than 2000 
steps. The implicit solution costs approximately 75% more 
per time step than the explicit solution because of the 
increased number of operations required to evaluate the flux 


jacobians. This increase in CPU time is far out-weighed by 
the faster convergence to steady state in that a converged 
solution is obtained with the implicit relaxation scheme with 
an order or magnitude less CPU time than the explicit scheme. 
The resulting steady pressure distribution is shown in Fig. 
3(b). For this case there is a relatively strong shock wave on 
the upper surface of the airfoil near 62% chord and a 
relatively weak shock wave on the lower surface near 30% 
chord. The pressure distributions indicate that there is only 
one grid point within the shock structure, on either the upper 
or lower surface of the airfoil, due to the sharp shock 
capturing ability of the upwind-type flow solver. 

Unsteady calculations with the upwind-type Euler solver 
were performed for the airfoil pitching harmonically about 

the quarter chord with an amplitude of a x = 2.51° and a 
reduced frequency based on semichord of k = 0.0814 at 

0.755 and a 0 ® 0.016°. These calculations are compared 
with the experimental data of Ref. 21. Instantaneous pressure 
distributions at eight points in time during the third cycle of 
motion using 2500 steps per cycle are shown in Fig. 4 for 
comparison with the experimental data. In each pressure plot 
the instantaneous pitch angle a(x) and the angular position in 
the cycle kx are noted. During the first part of the cycle there 
is a shock wave on the upper surface of the airfoil, and the 
flow over the lower surface is predominately subcritical. 
During the latter part of the cycle, the flow about the upper 
surface is subcritical, and a shock forms along the lower 
surface. The pressure distributions indicate that the shock 
position oscillates over approximately 25% of the chord along 
each surface, and in general, that the two sets of calculated 


Flux-vector splitting 


o Upper surface - Experiment 
o Lower surface - Experiment 





Fig. 4. Comparison of instantaneous pressure distributions for the NACA 0012 airfoil 
pitching atM_ = 0.755, a Q = 0.016°, = 2.51°, and k = 0.0814 computed 

using the upwind-type Euler flow solver with flux-vector splitting. 
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results compare well with the experimental data. Similar to 
the steady flow results, the shock waves are captured sharply 
with at most one grid point within the shock structure. The 
calculated results show the expected symmetry in the flow, in 
that the upper surface pressure distribution during the first 
half of the cycle is very similar to the lower surface pressure 
distribution during the second half of the cycle. A lack of 
similar symmetry in the experimental results suggests that 
the data may have been obtained at a slightly higher effective 
steady-state angle of attack than that reported in Ref. 21. 
Furthermore, the unstructured grid results of Fig. 4 are of 
comparable accuracy in comparison with published results 
obtained using structured grid methods for this case, such as 
those reported in Ref. 22. 


CflPCial-EuletfNavier-Stokes Results 

To demonstrate the central-difference-type conical 
Euler/Navier-Stokes flow solver and the adaptive mesh 
refinement procedures, calculations were performed for 
highly-swept delta wing and circular cone cases at high angle 
of attack and at supersonic freestream flow conditions. 13 In 
these calculations three cases were considered. In the first 
two cases, the conical Euler equations were solved for a 75° 
swept flat-plate delta wing at a freestream Mach number of 
1.4. In the first case, steady results were obtained at a - 20° 
angle of attack and p - 10° yaw angle. In the second case, 
unsteady results were obtained for the wing undergoing a 
forced harmonic rolling at a - 20° and P - 0°. For the third 



Fig. 5. Effects of adaptive mesh refinement on the total pressure loss contours for a 
75° swept flat-plate delta wing at - 1.4, a - 20°, and p - 10°; computed 
using the central-difference-type conical Euler flow solver. 
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case, the conical Navier-Stokes equations were solved for an 
85° swept circular cone at = 1.4, a = 20°, p - 0°, and 
Re = 0.5 x 106. 

Case_1. - The first case was selected to demonstrate 
applications of the adaptive mesh refinement procedures to a 
wing that is yawed. The calculations were performed by 
starting with a coarse mesh that has 32 nodes around the wing 
and 16 nodes in the outward direction for a total of 512 nodes. 
A partial view of this mesh is shown in the upper left part of 
Fig. 5. Results were obtained using adaptive mesh refinement 
by starting with the coarse mesh containing 512 nodes and 
adapting to the total pressure losses of the instantaneous 
solution to determine a locally fine mesh. The flow solver was 


run for a total of 4,000 iterations and the mesh was adapted at 
iterations 500, 1000, and 1500. Adapting the mesh every 
500 iterations produces intermediate results that are 
converged to plotting accuracy, although this is not a 
requirement for the adaption procedure to be numerically 
stable or to produce accurate final results. A summary of the 
meshes and the corresponding total pressure loss contours is 
shown in Fig. 5. For this case there is a strong flat vortex 
produced by the separated flow from the windward (left) 
leading edge with a crossflow shock wave beneath the vortex. 
There is also a weaker more-circular vortex produced by the 
separated flow from the leeward (right) leading edge. With 
the coarse grids containing 512 (original mesh) and 807 (one 
enrichment) nodes, the shock wave beneath the windward 



Fig. 6. Instantaneous total pressure loss contours during a cycle of harmonic rolling 
for a 75° swept flat-plate delta wing at - 1.4, a - 20°, p = 0°, o = 30°, 
and k = 0.3; computed using the central-difference-type conical Euler flow 
solver. 
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vortex is not detected in the total pressure loss contours. Upon 
further enrichment to a grid containing 1363 nodes, the 
shock begins to appear, and with the final grid of 2801 nodes 
the shock is shown to be sharply captured. The total pressure 
loss contours of Fig. 5 indicate that the vortical flow features 
are more clearly defined when the mesh is adaptively refined. 
A solution of comparable accuracy on a globally fine mesh 
would require 32768 nodes (256 x 128). A spatially 
accurate solution is thus obtained for this case with adaptive 
mesh refinement, using an order of magnitude fewer grid 
points (2801 nodes). 


Case 2. - The second case was selected to demonstrate the 
time-accurate capability of the conical Euler flow solver 
although the results are actually only locally conical. The 
unsteady results were obtained for the rolling delta wing 
oscillating harmonically with an amplitude of 30° at k - 0.3 
using 2400 steps per cycle of motion. In these calculations 
the grid was first adapted to the steady solution to create a fine 
embedded mesh locally which was then used in the unsteady 
calculation. The fine embedded region was made large enough 
so that the moving vortices and shock waves are always 
contained within this region and the total grid contained 5152 
nodes. It is recognized, however, that for unsteady 
applications a de-refinement procedure can be used in addition 



Fig. 7. Effects of adaptive mesh refinement on the total pressure loss contours for an 
85° swept circular cone at - 1.4, a - 20°, and Re - 0.5 x 10 6 ; computed 
using the central-difference-type conical Navier-Stokes flow solver. 
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to the refinement procedure, to remove as well as to add nodes 
and elements during the calculation. Also, during the unsteady 
calculation the grid was moved as a rigid body to conform to the 
instantaneous position of the wing. Three cycles of motion 
were computed to obtain a periodic solution. Instantaneous 
total pressure loss contours at eight points in time during the 
third cycle of motion are shown in Fig. 6. In each part of the 
figure, the instantaneous angular position in the cycle kt is 
noted where t is time normalized by one-half of the reference 
length and the streamwise freestream speed. 

During the first half of the cycle as the wing rolls 
clockwise to the maximum amplitude and then back to straight 
and level, the crossflow shock wave beneath the left vortex 
disappears, the left vortex weakens, and the right vortex 
grows in strength. During the second half of the cycle, the 
opposite situation occurs. The shock wave beneath the right 
vortex disappears, the right vortex weakens, and the left 
vortex grows in strength. It is easy to see that the solution is 
periodic in the third cycle of motion by comparing the total 
pressure loss contours at any two points in time that are 180° 
out of phase in the cycle. The two sets of total pressure loss 
contours, 180° out of phase, are antisymmetric as shown in 
Fig. 6. 



Case 3. - The third case was selected to further 
demonstrate the adaptive mesh refinement procedures for a 
viscous flow about a circular cone. The calculations were 
performed by starting with a coarse 64 x 32 mesh (2048 
nodes) and adaptively refining the mesh twice. The flow solver 
was run for a total of 4000 iterations and the mesh was 
adapted at iterations 1000 and 2000. A summary of the 
meshes and the corresponding total pressure loss contours is 
shown in Fig. 7. For this case, there are two large, 
symmetric, primary vortices that are produced by the 
separated flow from the left and right sides of the cone. With 
the coarse mesh containing 2048 nodes, the primary vortices 
are very diffuse, the secondary vortical flow features are 
absent, and the boundary layer appears thick. When the grid 


Experiment 
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Fig. 9. Comparison of lift coefficient versus angle of attack 
for the Langley fighter at M„ = 2.0 computed using 
the central-difference-type Euler flow solver. 


(a) original surface grid. 




a = 0° 

a « 12° 
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( b ) assumed bending mode. 

Fig. 8. Surface grid for the Langley supersonic fighter 
configuration. 


Fig. 10. Effects of angle of attack and reduced frequency on 
the lift coefficient responses due to the complete 
vehicle sinusoidal bending deformation of the Langley 
fighter configuration at M„- 2.0 computed using the 
central-difference-type Euler flow solver. 
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is refined to 3733 nodes and then to 8508 nodes, the primary 
vortices become smaller and more sharply defined, secondary 
and tertiary vortices become apparent in the flow, and the 
boundary layer appears to be much thinner in comparison 
with the coarse grid results. An accurate solution is thus 
obtained for this case with adaptive mesh refinement, using a 
factor of approximately four fewer grid points (8508 nodes) 
in comparison with a globally fine mesh (32768 nodes). 

Three-DimensionaLEuleLResults 

To assess the three-dimensional version of the central- 
difference-type Euler solver and dynamic mesh algorithm, 
calculations were performed for a supersonic fighter 
configuration that was tested at NASA Langley Research 
Center. 12 The fighter will hereafter be referred to as the 
Langley fighter 23 The Langley fighter is an aircraft designed 
to cruise at supersonic Mach numbers and also to manuever 
efficiently at transonic speeds. The aircraft has a high 
fineness-ratio fuselage with an underslung swept inlet and a 
flow-through duct. Other components which make up the 
configuration include a cranked wing with a leading edge sweep 
of 70° inboard and 20° outboard, a 55° swept canard, and a 
60° swept vertical tail. All of the lifting surfaces have thin 
circular-arc airfoil sections. Results were obtained for the 
Langley fighter using a grid which has 13,832 nodes and 
70,125 tetrahedra. The surface triangulation of the aircraft 
is shown in Fig. 8(a). Because of symmetry, the calculations 
were performed for only half of the aircraft- which contained 
4581 triangles on the surface of the vehicle. Steady-state 
results were obtained for the Langley fighter at a freestream 
Mach number of 2.0 and four angles of attack including a - 
0°, 4°, 8°, and 12°. Unsteady results were obtained at a « 0° 
and 12° for the aircraft oscillating harmonically in an 
assumed polynomial complete-vehicle bending mode which is 
shown in Fig. 8(b). Three values of reduced frequency based 
on wing tip semi-chord were considered including k= 0.025, 
0.05, and 0.1. The amplitude of the deformation was taken to 
be one-fifth of the mode shape shown in the figure. So, for 
example, the wing tip deflection in the calculation was 
approximately 20% of the tip chord. 

Steady flow results were obtained for the Langley fighter 
configuration using 1500 iterations starting from freestream 
flow conditions. Convergence to steady-state was accelerated 
using local time-stepping, implicit residua! smoothing, and 
enthalpy damping. In each case, the error in the solution, as 
defined by the L 2 -norm of the density residual, was reduced by 
approximately four orders of magnitude. Detailed color- 
contour plots of the pressure coefficients on the surface of the 
vehicle were presented in Ref. 12. A comparison between 



Fig. 11. Upper surface grid for the ONERA M6 wing. 


calculated and experimental lift coefficient versus angle of 
attack for the Langley fighter is shown in Fig. 9. In general, 
the calculated lift coefficient agrees well with the 
experimental data, especially at the lower angles of attack of 
0° and 4°. At the higher angles of attack of 8° and 12°, the 
calculated lift is slightly underpredicted in comparison with 
the measured values, which is possibly due to the coarseness of 
the grid. 

Unsteady flow results were obtained for the Langley fighter 
configuration oscillating harmonically in the complete-vehicle 
bending mode. The calculations were performed for k = 
0.025, 0.05, and 0.1 using 1440, 720, and 360 steps per 
cycle of motion, respectively, for five cycles. Calculated 
instantaneous pressure contours at the maximum (bend-up) 
and minimum (bend-down) amplitudes of oscillation, during 
the fifth cycle, were presented in Ref. 12. Lift coefficient 
responses for the fighter oscillating atM^- 2.0 are 
presented in Fig. 10. To quantify the effects of angle of attack 
and reduced frequency, five cycles of motion were considered 
as plotted. In each case, the mean value of the oscillating lift 
was subtracted off to allow for a direct comparison between 
results obtained at a - 0° and a » 12°. As shown in Fig. 10, 
the effects of angle of attack are very small for k = 0.025. As 
the reduced frequency is increased, however, the magnitudes of 
the response of the aircraft increase and there is a small 
change in phase of approximately 35°. Also, an effect due to 
angle of attack becomes evident, in that the responses at a = 
12° are of larger magnitude than those at a = 0°. This effect 
is possibly attributable to a small amount of vortical flow that 
is present along the upper surface of the wings in the inboard 
region at a = 12°. 



x/c 


Fig. 12. Comparisons of steady pressure distributions for the 
ONERA M6 wing at 0.84 and a Q « 3.06° 
computed using the upwind-type Euler flow solver 
with flux-vector splitting. 
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To test the more-recently-developed upwind-type Euler 
flow solver, calculations were performed for the ONERA M6 
wing. 24 The M6 wing has a leading edge sweep angle of 30°, 
an aspect ratio of 3.8, and a taper ratio of 0.562. The airfoil 
section of the wing is the ONERA "D" airfoil which is a 10% 
maximum thickness-to-chord ratio conventional section. The 
results were obtained using a grid which has 154,314 nodes 
and 869,056 tetrahedra. The surface triangulation for the 
upper surface of the wing is shown in Fig. 11. Results were 
obtained for the M6 wing at a freestream Mach number of 0.84 
and 3.06° angle of attack. These conditions were chosen for 
comparison with the experimental pressure data of Ref. 24. 
The results were obtained using the explicit time-marching 
scheme since it requires half of the memory of the implicit 
scheme. The code was run for 6000 time steps at a CFL 
number of 5.0, which produced a four order of magnitude 
reduction in the L. 2 -norm of the density residual. 

Figure 12 shows surface pressure coefficient comparisons 
with the experimental data at five span stations including T| - 
0.2, 0.44, 0.65, 0.9, and 0.95. In these plots the Euler 



(a) upper surface. 



(b) lower surface. 

Fig. 13. Surface pressure contour lines (Ap - 0.02) on the 
ONERA M6 wing at 0.84 and a Q « 3.06° 
computed using the upwind-type Euler flow solver 
with flux-vector splitting. 


results are given by the solid curves where plus signs have 
been included to indicate the actual grid point values which are 
connected with straight line segments. The experimental data 
is represented by the circles. For r| = 0.2 there are two shock 
waves along the chord. The forward shock wave is well 
predicted including the suction peak. The second shock wave is 
predicted slightly downstream of the experimental shock 
location which is typical of inviscid methods for this case. 
Also, the lower surface pressure coefficients agree well with 
the data. Att| - 0.44 the shock locations have begun to 
coalesce. The leading edge suction peak is well predicted and 
both shock waves are captured sharply. At q « 0.65 the 
forward shock wave is near 20% chord and the second shock 
wave is near midchord. All of the pressure levels are well 
predicted and both shocks are captured sharply with only one 
grid point within the shock structure. There are also no 
overshoots or undershoots near the shocks due to the flux 
limiting. Furthermore, the lower surface pressure 
coefficients are predicted accurately. At T] - 0.9 the two 
shocks have merged to form a single, relatively strong, shock 
wave near 25% chord. Here the shock is very sharply 
captured and the calculated pressures again agree well with the 
experimental data. Finally at T| = 0.95, the shock wave is 
slightly stronger than the previous span station. Here, the 
calculated shock again has only one interior point. 

Figure 13 shows pressure contour lines on the surface of 
the wing plotted using an increment of Ap = 0.02. Pressure 
contours on the upper surface are shown in Fig. 13(a); 
Pressure contours on the lower surface are shown in Fig. 
13(b). The upper surface contours (Fig. 13(a)) clearly show 
the lambda-type shock wave pattern formed by the two inboard 
shock waves which merge together near 80% semispan to form 
the single strong shock wave in the outboard region of the 
wing. The lower surface contours (Fig. 13(b)) indicate that 
there is very little spanwise variation in pressure. 


Concluding Remarks 

The current status of the development of unstructured grid 
methods in the Unsteady Aerodynamics Branch at NASA Langley 
Research Center was described. These methods are being 
developed for steady and unsteady aerodynamic applications. 
The paper highlighted the flow solvers that have been 
developed for the solution of the unsteady Euler and Navier- 
Stokes equations and gave selected results which demonstrated 
various features of the capability. The results demonstrated 
two- and three-dimensional applications for both steady and 
unsteady flows. Comparisons were also made with solutions 
obtained using a structured grid code and with experimental 
data to determine the accuracy of the unstructured grid 
methodology. These comparisons showed good agreement which 
thus verified the accuracy. 
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